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The  determination  of  gas  phase  thermal  rate  constants  for  chemical  reactions  even  at 
ambient  conditions  is  one  of  the  most  basic  problems  in  chemistry.  Thermal  rate  constants 
as  a  function  of  temperature  are  the  critical  component  of  all  atmospheric1  and  combustion 
models  24.  They  are  also  essential  to  understanding  the  initial  steps  in  detonation.  The 
enormous  number  of  potential  reactions— and  dearth  of  experimental  data  for  most  of 
them—  emphasize  the  critical  necessity  of  developing  procedures  that  can  reliably  provide 
a  priori  rate  constants  5  to  complement  experiment,  and  to  help  to  fill  in  the  gaps  in  the 
absence  of  experiment.  Despite  the  enormous  successes  of  quantum  mechanical 
descriptions  of  the  electrons  in  molecules--  i.e.  electronic  structure  theory-  where  several 
program  systems  exist  like  GAUSSIAN,  ACES  II  3,  MOLPRO  ,  TURBOMOLE,  QCHEM, 
GAMESS,  JAGUAR,  MOPAC,  HYPERCHEM  9,  etc,  that  are  extensively  used  to 
determine  the  structure,  spectra,  and  energetics  of  molecules-  the  equally  basic,  kinetic 
aspect  of  chemistry  has  not  been  as  developed.  Although  programs  and  program  suites 
like  POLYRATE  10,  CHEMRATE  ",  VariFlex  12,  UNIRATE  13,  and  MultiWell  1415  exist  for 
kinetics  calculations,  all  of  them  require  results  from  the  electronic  structure  theory  for 
input  information. 

As  we  have  pointed  out  in  our  proposal  entitled  “Computational  Prediction  of  Kinetic  Rate 
Constants”,  there  is  a  great  need  for  electronic  structure  methods  that  are  seamlessly  tied 
to  other  programs  that  permit  rate  constant  evaluation,  or  better,  have  accurate  rate 
constant  evaluation  as  an  integral  part  of  the  suite  of  electronic  structure  programs,  so  that 
rate  constants  can  be  evaluated  at  virtually  any  applicable  level  of  theory  by  non-expert 
users.  In  addition  to  ease  of  application,  further  attention  has  to  be  paid  to  improving  the 
electronic  structure  theory  that  enables  it  to  provide  the  highly  accurate  information 
specific  to  rate  constants,  that  their  reliable  evaluation  demands;  and  to  the  degree 
warranted,  improve  the  procedures  for  the  rate  constant  evaluation  itself. 

However,  the  challenges  of  a  priori  rate  constant  predictions  are  many: 

(1)  For  PES  with  intrinsic  barriers,  kinetic  theory  begins  with  transition  state  (TS) 
information,  defined  as  a  saddle  point  on  a  reaction  path  of  a  temperature 
independent  PES.  Unlike  molecules  at  normal  bond  lengths  and  angles,  TS’s  are 
characterized  by  beginning  to  break  and  form  bonds,  so  the  electronic  structure 
theory  description  is  more  demanding  in  terms  of  accuracy,  and  in  terms  of  the 
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such  as  activation  barriers  Once  this  is  done,  then  for  most  thermal  rate  constants  reasonable  values  can  be  obtained 
from  transition  state  theory  (TST),  variational  transition  state  theory  (VTST),  RRKM,  and  related  methods.  The  latter 
constitutes  the  dynamics  portion  of  the  effort,  but  the  computationally  determining  step  is  the  quantum  mechanical 
calculations  for  the  electronic  structure.  Most  of  the  accuracy  of  the  resultant  rate  constant  also  derives  from  that  part  of 
the  calculation.  Depending  upon  the  size  of  molecule,  there  is  a  transition  from  simpler  to  more  sophisticated 
methods. This  project  has  introduced  HF-DFT  as  the  simplest  possible  method  for  their  determination,  since  only  a  HF 
calculation  followed  by  a  ’tack-on’  DFT  estmate  of  correlation  is  required.  At  the  other  end  we  introduced 
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routine  application  There  is  also  an  interface  to  POLYRATE  to  further  exploit  the  ACES  II  input  to  that  package  for  rate 
constant  calculation.  Several  user  friendly  tools  were  added  to  ACES  II  to  make  such  rate  constant  applications  easier. 
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tools  (analytical  force  calculations  and  reaction  path  following)  needed  to  locate  the 
TS  s  on  the  PES. 

(2)  For  PES  without  intrinsic  barriers,  such  as  will  occur  with  recombination  reactions 
involving  two  radicals;  rather  than  the  energy  path  above,  it  is  necessary  to  find  the 
maximum  of  the  free  energy  path,  AG(T),  to  obtain  rate  constants  To  include  the  T 
dependence  the  appropriate  partition  functions  have  to  be  obtained  along  the  path, 
which  requires  vibrational  frequency  information  at  every  point. 

(3)  The  dynamics  part  of  the  problem  required  to  evaluate  the  rate  constants  once  the 
electronic  structure  problem  has  been  solved,  must  be  reliable  itself.  The  current 
choices  for  this  step  are  transition-state  theory  (TST),  RRKM,  variational  (VTST), 
semi-classical  generalizations  of  Miller’s  flux-flux  auto  correlation  function  approach 
(SCFFAF),16 17  classical  trajectories,  and  maybe  some  other  models  based  upon 
updated  collision  theory  concepts.  The  accuracy  of  the  results  obtained  will  depend 
greatly  upon  both  the  rate  evaluation  and  the  accuracy  of  the  electronic  structure 
results. 

(4)  Thermal  rate  constants  are  of  many  types  and  they  can  be  affected  by  tunneling 
when  light  atoms  are  involved,  van  der  Waals’  minima  on  the  PES18,  multiple 
reaction  channels  (transition  states),  and  non-adiabatic  effects  with  multiple  PESs. 
Tunneling,  in  particular,  is  usually  treated  as  an  ‘add-on’  to  VTST  with  various  small 
and  large  curvature  tunneling  corrections  possible.  The  FFAF  approach  7  begins  to 
put  tunneling  in  more  satisfactorily,  but,  in  practice,  only  within  a  semi-classical 
approximation  when  this  tool  has  been  used  for  multi-  atom  systems.18 

So  unlike  most  applications  of  electronic  structure  theory,  the  wealth  of  different  types  of 
reactions,  different  conditions,  and  the  more  global  PES  information  required,  makes  it 
extremely  difficult  to  develop  a  few,  general-purpose  procedures  that  can  be  refined  until 
accurate  rate  constants  are  obtained  from  something  approximating  a  ‘point  and  click’ 
computational  chemistry  calculation. 

During  the  three  years  covered  by  the  proposal  “Computational  Prediction  of  Kinetic  Rate 
Constants”  awarded  to  the  ACES  QC  and  the  University  of  Florida  sponsored  by  the 
United  States  Air  Force  Office  of  Scientific  Research  (AFOSR)  through  a  Phase  II  STTR, 
as  mandated,  we  have  primarily  focused  on  the  electronic  structure  development 
component  and  continued  developing  in  house  rate  program  (ACESRATE)  and  the 
interface  software  to  other  well  established  rate  programs  which  we  started  during  the 
phase  I  stage.  The  latter  developments  are  essential  in  order  to  make  the  electronic 
structure  data  easily  available  to  the  kinetic  programs.  Before  proceeding  to  discuss  them 
in  detail,  in  summary  our  accomplishments  include  the  following: 


Development,  implementation  and  calibration  of  two  new  electronic  structure 
methods  (FIF-DFT  and  ACCSD(T)). 
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•  Upgrading  KS-DFT  implementation,  adding  the  KS-DF  gradients  and  improving  the 
stationary  state  search  algorithms  in  ACES  II  and  in  ACES  III 

•  Developing  the  ACESRATE  program  suite  in-house  and  starting  to  develop 
interfaces  to  vender  supplied  rate  programs  such  as  POLYRATE  and  VENUS. 

•  Three  articles  in  peer  reviewed  journals  have  been  published  already  and  several 
more  are  in  preparation. 


Electronic  Structure  Method  Development 
ACCSD(T)  Energies  and  gradients: 

The  gold  standard  of  quantum  chemistry  is  CCSD(T)19  because  of  its  high  accuracy  near 
equilibrium  and  reasonable  computational  cost.  Unfortunately,  as  illustrated  in  Fig.  1 
below,  the  method  fails 


C arbon  M  onoxide 


dramatically  away  from  equilibrium  (using  a  spin-restricted  reference)  and  has  the 
incorrect  bond-breaking  slope  with  an  unrestricted  reference.  The  breakdown  of  CCSD(T) 
can  be  traced  to  its  dependence  on  CCSD  amplitudes.  For  a  bond-breaking  situation,  the 
CCSD  wavefunction  is  qualitatively  wrong;  however,  the  energy  is  only  a  linear  function  of 
the  amplitudes  and  is  well  behaved.  CCSD(T)  is  a  quadratic  function  of  the  CCSD 
amplitudes  leading  to  a  qualitatively  incorrect  energy. 

By  applying  similarity-transformed  perturbation  theory  to  the  CCSD  stationary  energy 
functional,  one  can  show  that  the  leading  order  correction  is  linear  in  the  T  amplitudes, 
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and  is  also  linear  in  the  A  amplitudes  from  CCSD  gradient  theory.  This  correction  is  called 
the  ACCSD(T)  method: 


E 


ACCSDtTl 


=  {O|0  +  A)H|O)+(||(«|AH 


ll)|  abr 
ijk 
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The  A  amplitudes  are  qualitatively  well  behaved,  and  that  substantially  improves  bond¬ 
breaking  behavior  for  the  perturbation  correction.  Only  twice  as  expensive  as  CCSD(T), 
ACCSD(T)20  represents  the  simplest  extension  (both  theoretically  and  computationally)  to 
CCSD(T)  that  improves  bond  breaking,  and  gives  results  of  the  same  quality  as  the  more 
complicated  and  formally  less  satisfactory  CR-CCSD(T)L  and  CCSD(2)T  methods.  The 
latter  are  either  not  size-extensive,  the  rationale  for  all  coupled-cluster  methods,  or  fail  to 
be  invariant  to  orbital  rotations. 


For  the  identification  of  stationary  points,  it  is  necessary  to  not  only  have  energies,  but 
also  analytical  derivatives.  For  a  perturbative  method  like  ACCSD(T)  this  is  not  trivial;  to 
be  able  to  apply  the  generalized  Flellman-Feynman  theorem,  a  stationary  formulation  of 
ACCSD(T)  needed  to  be  developed.  Then  the  derivative  is  defined  by, 

^£a«:si*T)  =(o|(|  +  A  +  n)hix\Q)  +  (0|(l  +  A^L^JO)  +  (0|A//"u'7f,|0) 

dx  +(0|Al{l/y<0,'r7^2||0)+(0|Al;l(^luT)r|0) 

We  have  now  implemented  ACCSD(T)  and  its  analytical  derivatives  to  make  such 
applications  routinely  possible  for  transition  states.  Having  correct  potential  energy 
surfaces  for  the  entire  range;  from  equilibrium  to  bond  breaking,  is  an  essential 
component  of  the  rate  constant  calculations.  We  found  that  despite  the  substantially 
increased  improvement  in  bond  breaking,  for  the  determination  of  stationary  points, 
including  transition  states,  CCSD(T)  and  ACCSD(T)  performed  virtually  identically21,22. 
The  transition  state  structure  of  the  nitromethane-methylnitrite  rearrangement  (in  Table  1) 
and  the  reaction  path  (Fig.  2)  computed  at  CCSD(T)  and  ACCSD(T)  support  this 
conclusion.  Another  system  that  we  studied  is  the  concerted  bond  breaking  of  RDX.  It  is 
shown  that  similar  to  nitromethane-methylnitrite  rearrangement,  both  the  CCSD(T)  and 
ACCSD(T)  show  essentially  the  same  characteristics.  As  shown  in  Table  2,  the  same 
conclusion  is  reached  for  the  barrier  heights  of  the  series  of  reactions  from  a  selection 
compiled  by  Truhlar  eta/.23 for  the  purpose  of  benchmarking  electronic  structure  methods 
for  rate  constant  calculations. 
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Fig  2.  NMT  to  MNT  reaction  path  computed  by  CCSD(T)  and  ACCSD(T). 
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Table  2:  Barrier  Heights  for  well-characterized  reactions  from  Minnesota  data  set. 
Energies  and  errors  from  best  experimental  estimates  are  in  kcal/mol. 
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Non  Variational  Density  Functional  Theory  (HF-DFT)  Energies  and  Gradients: 


The  non  variational  DFT  (HF-DFT)24  methods  do  not  involve  solving  Kohn-Sham 
equations.  They  are  designed  to  replace  the  HF  exact  exchange  energy  by  the  energy  of 
a  chosen  exchange-correlation  functional  computed  with  the  HF  density.  The  assumption 
is  that  by  choosing  an  exchange  correlation  functional  instead  of  the  HF  exchange,  the  HF 
energy  is  improved  by  the  correlation  effects  introduced  by  the  exchange  correlation 
functional.  The  HF-DFT  energy  and  gradient  expressions  are  straightforward  extensions 
the  HF  energy  and  gradient  expressions.  For  example,  the  energy  and  the  gradient 

dFHF~DFT 
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(<P  )  +  ££( 


c-m* 
change  correlation 


(P 


liar  tree-  hock 


:kyp 


Hartree-Eock 
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/W//-/WT 

db  _  yy  pus  a  ai  v< 


dR, 


(7  / IV 


dR 


+  - 


d£ 


III  l)IT 


dR 


pO  _  "SP  p'CT/^CT 

Where  the  density  matrix  jjv  Z-*  /"  m  and  the  KS-Fock  operator  (for  the  gradients) 

teocc 

is  given  by 


F, 


AY  <7 

HV 


I  M  7 


XS 


When  expanded  the  HF-DFT  energy  gradient  expression 

t^ZVC  +t  +£:APa,pp) 


dx 


1  *— 

pvAo 


-ZZZ^^ZZZ":^ 

(T  t  )  (T  i  a 

Where  En  =\ f(ptrpe . )dr  ^  =SZ wsuf(P? . *4) 

i 

^  II  Ex,  =  Z  Z  l  U’.W  V  $f(Pa< '/W  )  +  V II  W.M  f(p“ r,\i)  ] 


For  the  most  part,  the  cost  involving  adding  the  exchange  correlation  energy  is  small  and 
one  can  safely  assume  that  there  is  no  additional  cost  beyond  the  cost  of  obtaining  the  HF 
solutions,  but  the  energies  are  corrected  for  correlation.  Furthermore,  since  it  depends 
upon  a  HF  solution  to  get  the  density,  there  is  no  self-interaction  failure,  but  all  correlation 
effects  are  assumed  to  come  solely  from  the  evaluation  of  the  LYP  functional,  in  a  non¬ 
self-consistent  manner.  The  results  for  a  large  number  of  molecules  have  been 
surprisingly  accurate,  as  the  mean  absolute  errors  (MAE)  compare  favorably  to  those  from 
the  very  high-level  CCSD(T)  method.  For  example  for  equilibrium  geometries,  the  average 
absolute  error  for  a  large  number  of  small  molecules  compared  to  experiment  is  0.005 
Angstroms  for  bond  lengths,  the  same  as  CCSD(T),  1 .7  degrees  for  the  bond  angle, 
compared  to  1 .9  for  CCSD(T),  40  cm  error  in  vibrational  frequencies,  compared  to  30  cm 
1  for  the  latter,  and  does  better  for  atomization  energies,  3.6  kcal/mol  compared  to  1 1 .5 
kcal/mol24. 


The  primary  data  that  are  of  interest  for  rate  constant  calculation  are  activation  energies, 
structures,  and  vibrational  frequencies  of  stationary  points  (and  in  some  cases  at  points  in 
the  vicinity  of  a  transition  state).  Until  this  work,  the  performance  of  HF-DFT  for  transition 
state  structures,  frequencies  and  activation  barriers  was  not  known.  To  facilitate  such 
calibrations,  we  implemented  the  analytical  gradients  for  HF-DFT^.  In  Table  3,  we  show 
the  structures  of  the  transition  states  of  a  series  of  reactions  (18  altogether)  obtained  with 
HF-DFT  along  with  the  QCISD  results  from  the  literature  obtained  with  the  same  MG326 
basis  sets.  As  shown  in  Table  4,  the  average  absolute  difference  between  the 
QCISD/MG3  and  HF-DFT/MG3  results  are  0.03  (Angstroms)  and  2.6  (degrees)  for  bond 
lengths  and  angles  respectively. 
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The  accuracy  of  the  activation  energy  is  perhaps  the  most  important  factor  that  governs 
the  accuracy  of  the  thermal  rate  constants  since  the  error  in  the  rate  constant  grows 
exponentially  with  the  error  in  the  activation  energy.  In  order  to  calibrate  the  accuracy  of 
HF-DFT  activation  energies,  we  have  done  two  sets  of  calculations.  In  one  case,  we  have 
used  the  geometries  optimized  at  the  MG3/QCISD  level  available  in  the  literature  and  the 
single  point  HF-DFT  energies  of  reactants  and  transition  states  are  computed.  This  set  is 
comprised  of  40  reactions  as  compiled  by  Truhlar  and  coworkers23  in  order  to  calibrate 
new  DFT  functionals.  In  the  other  set  (a  subset  of  reactions  from  the  previous  set)  the 
activation  energies  are  computed  at  the  HF-DFT  optimized  transition  state  structures.  For 
comparisons,  we  also  show  the  HF  and  KS-DFT  activation  energies  (also  obtained  at  the 
TS  geometry  optimized  at  the  HF  and  KS-DFT  level).  As  shown  in  Table  5,  the  average 
absolute  error  (the  best  available  values  from  the  literature  is  chosen  as  the  reference)  of 
the  HF-DFT  activation  energies  is  2  kcal/mol.  A  2  kcal/mol  average  error  for  HF-DFT  is 
remarkable  in  the  sense  that  there  is  practically  no  extra  cost  beyond  the  cost  for  HF 
calculation.  That  is,  this  is  the  simplest  possible  calculation  that  includes  electron 
correlation,  a  necessity  for  activation  barriers  and  transition  states.  We  conclude  that  HF- 
DFT  is  a  cost  effective  and  sufficiently  accurate  method  to  be  used  in  rate  constant 
calculations. 


Table  3:  Comparisons  of  HF-DFT  transition  state  structures  with  publish  data  (distances 
are  in  angstroms  and  angles  are  in  degrees). 


reaction  QCISD/MG3 

HFDFT/MG3 

error 

Abs(error) 

OH  +  H2  ->  H  +  H20 

t.l;.  -  .  x. . 

tst(H-O.-.H-H) 

t 

R(H-O) 

0.967 

0.978 

-0.01 

0.01 

R(O...H) 

1.3165 

1.4121 

-0.1 

0.1 

R(H-H) 

a(h~o...h) 

0.833 

0.805 

0.03 

0.03 

97.19 

97.96 

-0.77 

0.77 

a(6...h-h) 

165.48 

160.46 

5.02  5.02 

CH3  +  H2  ->  H  +  CH4 

tst(CH3....H-H) 

_ _ _  _  ...  1  _ 

R(C...H) 

1.382  1.431 

-0.05 

0.05 

R(C-H1) 

1 .083  1 .086 

-0  0 

■  _  .  ..  __  _ _ _ _  ... 

R(H-H) 

0.9040  0.8773 

0.03  0.03 

...  -  t —  -  - .  . 

A(C..H-H) 

180  180 

0  0 

A(C-H..H) 

114.75 

102.49 

12.3 

12.3 

OH+CH4  ->  CH3  +  H20 

tst(CH4..0-H) 

R(O-H) 

0.966 

0.977 

-0.01 

0.01 

R(H...O) 

1.227 

1.396 

-0.17 

0.17 
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R(C-H)  1 .084 

1.088 

-0 

0 

A(C-H..O)  173.68  174.08 

-0.4 

0.4 

A(H....O--H)  97.7  97.38 

0.32 

0.32 

A(H...C--H) 

104.25  104.65 

-0.4 

0.4 

A(H-C-H) 

112.99  112.33 

0.66 

0.66 

H  +  CH30H  ->  CH20H  +  H2 

tst(HOCH3..H) 

R(c--O) 

1.384 

1 .396 

-0.01  I 

0.01 

R(O--H) 

0.957 

0.967 

-0.01 

0.01 

R(C--H) 

1.089  1.095 

-0.01 

0.01 

R(C-H) 

1.084  1.089 

-0.01 

0.01 

R(C...H) 

1.315  1.294 

0.02 

0.02 

R(H....H) 

0.969  0.987 

-0.02 

0.02 

A(C--O-H) 

109.207 

109.78 

-0.57 

0.57 

A(O-C--H) 

115.07 

115.04 

0.03 

0.03 

A(C....H..H) 

177.479 

178.16 

-0.68 

0.68 

A(H-CH) 

109.207 

109.78 

-0.57 

0.57 

A(H-CH) 

109.819 

109.58 

0.24 

0.24 

A(H-CH) 

110.068 

110.65 

-0.58 

0.58 

OH  +  NH3  -  >  H20  + 
NH2 

tst(NH3..0H) 

1  ’ 

R(N--H) 

1.016  1.020  -0 

0 

R(N-H) 

1.016  1.020  -0 

0 

R(N-H) 

1.153  1.100 

0.05 

0.05 

R(H....O) 

1 .259  1 .383 

-0.12 

0.12 

R(H-O) 

0.964  0.975 

-0.01 

0.01 

A(H-N-H) 

106.64  107.21 

-0.56 

0.56 

A(H-N-H) 

108.25 

106.73 

1.53 

1.53 

A(N-H-O) 

150.42 

148.05 

2.37  2.37 

A(H--O-H) 

1 - - - *— **  '  1  *'  f  “'"rir“r'  ' 

104.31 

106.8 

-2.48  2.48 

F  +  H2  ->  H  +  HF 

 I 

tSt(H--H.-.F) 

R(H--H) 

0.774 

0.7707 

i  -  ..  . _ ,  ......  .  . 

0 

0 

R(H...F) 

1.466 

1 .5886 

-0.12 

0.12 

A(H-H...F) 

152.06 

108.16  43.9 

43.9 

H  +  PH3  ->  PH2  +  H2 

i _ _ _ . _ . _ _ „ _ 

tst(PH2~H...H) 
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H  +  H2S  ->  H2  +  HS 


R(H-S) 

1.337 

1 .3487 

R(S-..H) 

1.429 

1.3965 

R(H...H) 

1.149 

1.332 

A(H-S-..H) 

91.152 

91.884 

A(S-..H...H) 

172.61 

179.3 

-0.01 


0.03 


-0.18 


-6.68 


R(P-H) 

1.4139  1.424  -0.01 

0.01 

R(P-H) 

1.414  1.424  -0.01 

0.01 

R(P-..H)  1.484  1.46  0.02 

0.02 

R(H....H)  1.264  1.475  -0.21  0.21 

A(H-P-H) 

93.49 

92.908  0.59 

0.59 

A(H-P-..H) 

92.57  92.397 

0.17 

0.17 

A(P-..H..H) 

166.75  171.12 

-4.37 

4.37 

H  +  CIH'  ->  HCI  +  H 

tst(H-CI....H) 

R(H-CI) 

1 .485  1 .477 

0.01 

0.01 

R(CI...H) 

1 .485  1 .478 

0.01 

0.01 

A(H--CL.H) 

180  180 

0 

0 

OH  +  H  ->  H2  +  0 

1  “  j  “ 

tst(0-H...H) 

R(O--H) 

1.189 

1.3128 

-0.12 

0.12 

R(H...H) 

0.911  0.8358 

1 -  -  -  -  -  -  L --4 — — 

0.08 

0.08 

A(0-H...H)  180  180 

0 

0 

f  .  —  ■  -  - 

H  +  trans-N2H2  ->  H2  +  N2H 

tst(HNN~H...H) 

R(N-N)  1.225 

1 .2582 

-0.03 

0.03 

R(H-N)  1.030  1.0185 

0.01 

0.01 

R(N-..H)  1.117 

1.0619 

0.06 

0.06 

R(H-H)  1.179 

1 .3798 

-0.2 

0.2 

A(H-N-N)  107.41 

121.81 

-14.4 

14.4 

A(N-H-H)  170.55  169.49 

1.07 

1.07 

A(N-N-..H)  107.04  119.43 

-12.4 

12.4 

0.01 


0.03 


0.18 


0.73 


6.68 
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0  +  HCI  ->  OH  +  Cl 

tst(O...H~CI) 

R(O...H) 

2.482  2.4829 

1.414  1.4706 

0 

0 

R(H-...CI) 

-0.06 

0.06 

A(O...H-..CI) 

22.43  19.435 

3 

3 

CH4  +  NH  ->  NH2  +  CH3 

tst(H3C-...H..N-H) 

R(H-C)  1.083  1.0867 

-0 

0 

R(C-..H)  1.3953 

1.4002 

-0 

0 

R(H...N) 

1.2117 

1.2059 

0.01 

0.01 

R(N-H) 

1.0315 

114.59 

1 .0387 

-0.01 

0.01 

A(H-C-H) 

114.72 

-0.12 

0.12 

A(H-C-..H) 

101.54 

114.72 

-13.2 

13.2 

A(H...N-H) 

174.46 

176.36 

-1.9 

1.9 

- •  - -  I 

C2H6  +  NH  ->  NH2  + 

C2H5 

tst(H3C-CH2-..H...N-H) 

R(H-C) 

1.093 

1 .0987 

-0.01 

0.01 

R(H-C) 

1.090 

1.094 

-0 

0 

R(H-C) 

1.089 

1.094 

-0 

0 

R(C-C) 

1.508 

1.5141 

-0.01 

0.01 

R(C-H) 

1.086 

1.0912 

-0 

0 

R(C-H) 

1.086 

1.0911 

-0 

0 

R(C-..H) 

1.372 

1.337 

0.04 

0.04 

R(H...N) 

1.233 

1 .2639 

-0.03 

0.03 

R(N-H) 

1.032 

1 .0386 

-0.01  0.01 

A(H-C-H) 

107.447 

107.3 

0.14  0.14 

A(H-C-H) 

107.451 

107.86 

-0.41 

0.41 

A(H-C..C) 

1 1 1 .006 

107.86 

3.15 

3.15 

A(C--C-H) 

115.596 

115.55 

0.04 

0.04 

A(C-C-H) 

115.75 

115.54 

0.21 

0.21 

A(C-C-..H) 

106.615 

105.78 

0.84 

0.84 

A(C-  .H..N) 

175.085 

178.57 

-3.49 

3.49 

A(H...N-H) 

94.92 

99.646 

-4.73 

4.73 

-  ...  . — 

C2H6  +  NH2  ->  C2H5  +  NH3 

tst(H3C-CH2- 

..H..NH2) 

R(H-C)  1.093  1.0975 

-0 

0 
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Cl  +  H2  ->  HCI  +  H 

QCISD 

Best- 

QCISD 

HFDFT 

i . - . .  , 

Best- 

HFDFT 

QCISD- 

HFDFT 

Best 

tst(CI...H — H) 

.  i   ... . 

R(CI...H) 

1.543 

-0.113 

1.361 

0.069  0.18 

1.43 

R(H--H) 

0.774 

0.207 

1.184 

-0.203  -0.41 

0.981 

*The  worst  example,  presented  first  in  Table  1 ,  would  appear  to  be  Cl  +  H2,  but  this  is 
misleading.  In  this  case  the  QCISD/MG3  ‘reference’  structure  is  very  poor.  So  actually,  the 
error  for  HF-DFT  compared  to  the  best  possible  results  is  much  better,  attesting  to  HF- 
DFT’s  accuracy  (the  “best  values”  are  taken  from.  Phys.  Chem.  A  2001,  105,  2936-2941) 


Table  4:  The  error  analysis  of  HF-DFT  bond  lengths,  angles  and  activation  energies  (the 
reference  values  are  taken  from  J.  Phys.  Chem.  A,  Vol.  104,  4811, 2000).  The  distances, 
angles  and  energies  are  in  angstroms,  degrees  and  kcal/mol  respectively.) 


Error 

(QCISD-HF-DFT) 

ABS(QCISD-HFDFT) 

Distance(angstroms) 

-0.02 

0.038 

Angle(degree) 

0.1 

2.64 

Errors 

Best-HFDFT 

Best-QCISD 

HFDFT 

QCISD 

- 

(methodl  -method2) 

-0.12 

-1.02 

-2.13 

Absolute(method1- 

method2) 

2.28 

2.47 

3.17 
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Table  5:  Comparison  of  the  HF-DFT  activation  energies  with  the  literature  values.  The  HF- 
DFT  activation  energies  are  computed  at  the  same  geometry  (QCISD/MG3)  using  the 
same  basis  set  (MG3)  and  the  literature  data  is  complied  by  Truhlar  and  coworkers23. 


reactions  HFDFT/MG3//QCISDM33  QC 1$  D/MG3  //  QCI5D/TVIG3  best  activation  energy  HF  DFT*QC  IS  C  HF  DFT  Bestvalue 

kcal/mol  kcalAnol  kcalAnol  kcalAnol  kcalAnol 


CI+  H2  ->  HCI+  H 

S.92 

12  73 

8.7 

6.81 

2.78 

HfHCl  ->  Cl  +  H 2 

1.24 

7.01 

5.6 

-5.77 

-4  36 

OH  +  H2  ->  H  +  H20 

4  91 

7.77 

5.7 

-2  86 

-0.79 

H  +  H20  >  OH  +  H2 

18.42 

21.2 

22 

2  78 

-3  58 

CH3  +  H2  ->  H  +  CH4 

10.64 

13.82 

12.1 

-3.18 

-1.46 

CH4  +  H  ->  CH3  +  H2 

10.97 

15.75 

15 

-4.78 

-4.03 

OH  +  CH4  ->  CH3  +  H20 

S  87 

9.34 

6.7 

-3.47 

-0.83 

H20  +  CH3  ->  OH  +  CH4 

19  06 

20,84 

20  2 

-1.78 

-1.14 

H  +  CH30H  ->  C  H20H  +  Hi 

5. 58 

10.96 

7.3 

-5  38 

-1.72 

H2  +  CH20H  ->  H  +  CH30H 

14  63 

17.49 

13  8 

-2.86 

0  83 

H  +  h2  ->  H2  +  H 

5.23 

10.68 

9.6 

-5.45 

-4.37 

OH  +  NH3  -  >  H20  +  NH2 

3.3 

7.55 

3.2 

-4.25 

0.1 

NH2  +  H20  ->  NH3  +  OH 

13  7 

17.07 

13.2 

-3.37 

0.5 

HCI+  CH3  >  CI+  CH4 

1.08 

3.72 

1.8 

2.64 

-0.72 

CH4  +  Cl  >  CH3  +  HCI 

6.09 

11.37 

7.8 

-5.28 

-1.71 

OH  +  C2H6  ->  H20  +  C2HS 

4.4 

646 

3.4 

*2.06 

1 

H20  +  C2HS  >  OH  +  C2H6 

21  89 

2122 

20.7 

0  67 

1.19 

IF  +  H2  ->  H  +  HF 

1.83 

3.31 

1.8 

-1.48 

0.03 

H  +  HF  ->  F  +  H2 

32  95 

31  95 

33.2 

1 

-0  25 

OH  +  CH3->  0  +  CH4 

8  26 

10.32 

78 

-2.06 

0.46 

0  +  CH4  ->  OH  +  CH3 

9.37 

17.84 

13.7 

-8.47 

-4.33 

H  +  PH3  >  PH2  +  H2 

0.93 

4  09 

3.2 

-3.16 

-2.27 

H2  +  PH2  ->  H  +  PH3 

25.55 

27.43 

25.5 

-1,88 

0.05 

H  +  C IH  ->  HCI  +  H 

14.47 

20  66 

18 

-6.19 

-3  53 

OH  +  H  ->  H2  +  0 

7.52 

10.94 

10.1 

-3  42 

2  58 

0  +  H2  >  OH  +  H 

8.31 

16.53 

13.1 

-8.22 

-4.79 

H  +  trans -N2H2  >  H2  +  N2H 

8.24 

4.27 

5.9 

3.97 

2.34 

H2  +  N2H  >  H  +  trans  N2H2 

42.03 

42.88 

41.1 

-0.85 

0.93 

H  +  H25  ->  H2  +  HS 

1.21 

5.09 

3.6 

-3.88 

-2  39 

HS  +  H2  ->  H  +  H2S 

17.6 

21.15 

17.4 

-3.55 

0.2 

0  +  HCI  ->  OH  +  Cl 

9  92 

14.2 

9.8 

-4.28 

0.12 

Cl  +  OH  >  0  +  HCI 

13.81 

14  32 

9.9 

-0.51 

3.91 

CH4  +  NH  ->  NH2  +  CH3 

19.29 

25.39 

8.4 

6,1 

10.89 

CH3  +  NH2  ->  NH  +  CH4 

9.4 

10.67 

22.7 

-1.27 

-13.3 

C2H6  +  NH  ->  NH2  +  C2H5 

16.99 

22.46 

8 

-5.47 

8.99 

C2HS  +  NH2  ->  NH  +  C2H6 

11.4 

11 

18.4 

0.4 

-7 

C2H6  +  NH2  >  C2H5  +  NH3 

11.56 

14.48 

10  4 

-2  92 

1.16 

NH3  +  C  2H5  ->  C2H6+  NH2 

18.65 

19.71 

17.8 

-1.06 

0.85 

NH2  +  CH4  >  CH3  +  NH3 

13.5 

17.09 

14.5 

-3  59 

-1 

NH3  +  CH3  >  CH4  +  NH2 

16.29 

19  06 

17.9 

-2.77 

1.61 

CsH3^  C5H8 

34.76 

44.27 

38.4 

-9.51 

-3.64 
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Table  6:  A  comparison  of  HF-DFT  activation  energies  (in  kcal/mol)  with  HF  and  KS-DFT. 
The  structures  are  optimized  at  the  HF-DFT  level  using  the  MG3  basis  set. 


Reactions 

HF/MG3 

HF/MG3 

// 

HFDFT/MG3// 

HFDFT/MG3 

DFT/MG3//DF 

T/MG3 

Best 

Forward 

Back 

Forward 

Back 

Forward 

Back 

Forward 

Back 

CH3  +  H  2 

H  +  CH4 

22.1 

23.6 

10.8 

11.1 

7.2 

7.8 

12.1 

15 

H  CH3OH 

ch2oh  +  h2 

20.1 

24.2 

5.7 

14.7 

1.0 

11.7 

7.3 

13.8 

F  +H2 

H  +  HF 

-0.1 

14.1 

2.0 

33.2 

0.0 

32.3 

1.8 

33.2 

H  +PH3 

PH2  +  h2 

11.5 

31.4 

1.7 

26.2 

0.0 

24.2 

3.2 

25.5 

H  +CIH* 

HCI  +  H' 

31.5 

31.5 

14.5 

14.5 

9.8 

9.8 

18 

18 

OH  +  H 

h2  +  o 

18 

32.9 

9.4 

10.1 

1.6 

1.0 

10.1 

13.1 

H  +  H2S 

H2  +  HS 

12.7 

27.8 

2.0 

18.4 

0.0 

16.1 

3.6 

17.4 

0  +  HCI 

OH  +  CI 

28.3 

20 

11.8 

15.7 

-3.2 

2.0 

9.8 

9.9 

CH4  +  NH 

nh2+  ch3 

39.3 

20.5 

19.3 

9.5 

13.6 

3.6 

22.7 

8.4 

C2H6  +  NH 

nh2  +c2h5 

37.4 

21.1 

16.8 

11.2 

10.1 

5.2 

18.4 

8 

C2H6  +  NH. 

2  nh3+c2h5 

30.1 

30.7 

12.1 

19.1 

5.2 

12.8 

10.4 

17.8 

NH2  +  ch4 

ch3  +  nh3 

31.9 

30 

13.7 

16.5 

7.8 

10.4 

14.5 

17.9 

Table  6-1 :  Error  analysis  of  the  activation  energies. 


Mean  Errors 

(Method-Best) 

Kcal/mol 

Absolute(Method-Best)  Kcal/mol 

Hartree-Fock 

12.6 

12.9 

HF-DFT 

-1.0 

1.6 

DFT 

-6.6 

6.6 

Table  6-2: .  Error  analysis  of  above  12  reaction’s  Transition  state  geometries. 


Mean  Errors 

(QCISD-method)  Bond 
distance  in  Angstrom 

(QCISD-method)  Bond  angle  in 
degree 

Hartree-Fock 

0.1 

5.2 

HF-DFT 

0.02 

2.6 

DFT 

0.11 

0.5 
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Enhancing  the  Capabilities  of  ACES  II 

Development  and  Implementation  of  the  Intrinsic  Reaction  Coordinate  (IRC): 

We  have  implemented  the  state-of-the-art  IRC  method  into  ACESRATE  to  map  the 
reaction  paths.  In  this  module,  the  IRC  is  calculated  using  the  quantum  mechanical  data 
generated  by  means  of  ACES  II  and  ACES  III.  The  commonly  used  two  algorithms  have 
been  implemented: 

(1)  Euler  algorithm, 

(2)  Page-Mclver  algorithm, 

as  well  as  the  newly  developed  correction  method  on  top  of  those  two  kinds  of  IRC  paths: 

(3)  Euler  and  Page-Mclver  algorithm  with  correction. 

The  Euler  method  (1)  needs  only  the  energy  and  gradient  from  quantum  chemical 
calculations  in  each  step  and  thus,  is  feasible  without  heavy  quantum  chemical 
calculations.  However,  to  obtain  the  highly  accurate  IRC,  which  is  mandatory  for  the  rate 
constant  calculations,  a  very  small  step  size  is  required  to  accurately  follow  the  IRC 
causing  a  prohibitively  large  number  of  quantum  chemical  calculations  to  be  necessary. 
The  Page-Mclver  method  (2)  improves  the  accuracy  of  following  the  IRC  by  requiring,  in 
addition,  the  information  from  the  Hessian.  Consequently,  the  Page-Mclver  method  is 
able  to  use  a  larger  step  size  with  maintaining  the  accuracy  of  the  IRC  and  to  reduce  the 
number  of  points  on  which  quantum  chemical  calculations  are  performed.  Although,  the 
calculation  of  the  Hessian  itself  is  time  consuming,  the  Page-Mclver  method  should  be 
preferred  since  for  the  variational  transition  state  theory  (VTST)  rate  constant  calculations, 
the  calculation  of  the  Hessian  is  compulsory  anyway. 

In  the  course  of  several  application  calculations,  we  have  found  typical  examples  in  which 
the  Page-Mclver  method  produces  an  inadequate  IRC.  Figure  1  shows  the  modeled 
potential  surface  of  the  dissociation  reaction, 

CH3POF2  -»  CH3  +  POF2,  (1) 

in  which  the  reactant  has  a  staggered  structure  while  the  products  favor  the  eclipsed 
conformation.  IRC  is,  theoretically,  the  path  which  connects  the  reactants  and  products 
smoothly  while  keeping  the  correct  Hessian  structure  (all  the  eigenvalues  which 
correspond  to  the  modes  normal  to  the  IRC  are  positive)  as  shown  by  the  dashed  line  in 
figure  1.  The  solid  line  is  the  path  predicted  by  the  Page-Mclver  method.  Along  this  path, 
the  Hessian  does  not  hold  its  correct  structure,  which  makes  the  calculation  of  the  rate 
constant  unfeasible.  We  have  found  that  this  is  because  the  gradient  component 
corresponding  to  internal  rotation  is  very  small  relative  to  the  component  corresponding  to 
dissociation  and  Page-Mclver  method  tends  to  ignore  the  small  component.  We  have 
modified  the  method  to  correct  the  direction  of  the  Page-Mclver  step  by  partitioning  the 
Hessian  and  gradient  into  the  IRC  direction  and  its  orthogonal  counterparts  and  by 
applying  energy  minimization  within  the  orthogonal  space  to  take  full  account  of  the  small 
gradient  component. 

Figure  3  shows  the  IRC’s  produced  by  the  Page-Mclver  method  (red  line)  and  by  our 
correction  (blue  line.)  It  is  obvious  that  the  improvement  of  IRC  by  the  new  correction 
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method  is  remarkable  when  the  Page-Mclver  method  is  incapable  of  generating  an  IRC 
even  with  a  very  small  step  size  (0.005  A  amu  1 2.)  Figure  4  also  shows  the  effect  of  step 
size  used  for  the  Page-Mclver  with  correction  method.  The  step  size  of  0.1  A  amu  1 2  is 
the  choice  which  is  commonly  accepted  for  most  of  the  latest  IRC  algorithms.  As  shown 
in  figure  2,  all  the  step  size  tested  in  this  study  produce  remarkably  accurate  IRC’s  up  to  a 
step  size  of  0.5  A  amu  1/2,  which  is  extraordinary  large.  Table  I  collects  the  rate  constants 
estimated  from  the  IRC's  calculated  with  the  respective  step  sizes.  It  shows  that  all  the 
rate  constants  agree  very  well  each  other  regardless  of  the  step  size  for  the  respective 
temperature.  Since  the  larger  step  size  necessitates  fewer  quantum  chemical 
calculations,  it  is  evident  that  the  new  correction  method  enables  the  accurate  yet  efficient 
calculations  of  the  IRC  which  is  vital  for  accurate  yet  convenient  rate  constant 
calculations. 
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Figure  3.  Schematic  representation  of  the  potential  Figure  4.  IRCs  calculated  by  means 
surface  for  the  reaction  (1).  The  solid  line  is  the  IRC  °f  the  Page-Mclver  method  (red)  and 
calculated  with  the  Page-Mclver  method  while  the  our  correction  (blue), 
dashed  line  is  the  correct  IRC. 


Integration  of  KS-DFT  to  ACES  II  and  Implementation  of  KS-DFT  Analytic  Gradients: 

Not  only  were  the  DFT  capabilities  that  were  in  existence  in  ACES  II  underperforming 
but  they  also  they  were  not  properly  integrated  to  the  rest  of  the  ACES  II  program  system. 
These  two  limitations  were  severely  hampering  the  further  development  of  DFT 
capabilities  and  implementing  gradients,  Hessians  and  higher  derivatives  on  one  hand, 
and  on  the  other,  limiting  the  reactions  that  could  be  investigated  to  small  molecules.  In 
order  to  remedy  these  limitations  we  have  completed  a  software  integration  and 
improvement  of  the  DFT  programs  in  ACES  II  system.  As  a  result  of  these  improvements, 
the  current  DFT  calculations  take  2-3  times  less  than  the  previous  version  and  can  be 
distributed  as  part  of  ACES  II.  Also,  as  we  have  shown  in  Table  7,  that  benchmark 
calculations  of  B2H6  show  that  ACES  II  DFT  performance  is  within  a  factor  of  three  on  a 
per  iteration  basis,  from  the  NWChem  program  system.  The  time  taken  per  iteration  is 
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given  in  square  brackets  with  the  number  of  iterations  in  parentheses.  It  is  widely  thought 
that  for  the  standard  DFT  methods  NWChem  provides  the  fastest  implementation.  As  our 
DFT  methods  include  some  that  are  not  standard,  FIF-DFT  and  ab  initio  dft,  small  timings 
differences  will  likely  be  less  significant  to  our  objectives. 


Table  7:  The  performance  comparison  of  ACES  II  and  NWChem  B2H6  KS-DFT 
implementations  (without  symmetry;  Sun  OS). 


Number  of  functions 
and  memory 

cc-pVDZ  58,15MB 

cc-pVTZ  144,15MB 

aug-cc-pVTZ  230,15MB 


ACESII  NWChem 

Time  (h:m:s) 


0:2:38  (22)  [5] 
0:22:25  (23)  [58] 

1 :27:40  (23)  [228] 


0:0:24(8)  [3] 
0:1:46(8)  [13] 
0:10:26(8)[78] 


B2H6  (with  symmetry;  Sun  OS) 

aug-cc-pVTZ  230,15MB  _  0:28:34(21 )  [81]  0:2:34(2)  [77] 


The  analytic  gradient  capabilities  for  KS-DFT  are  also  implemented.  In  Table  8,  we  show 
preliminary  timing  data  for  ACES  II  and  NWChem  KS-DFT  single  point  gradient 
calculations.  Both  C6H6  and  B2H6  use  D2h  symmetry  and  in  both  cases,  NWChem 
performs  better  than  ACES  II.  However,  for  the  N10  molecule  which  uses  no  symmetry, 
ACES  II  outperforms  NWChem.  While  both  ACES  II  and  NWChem  carry  out  the  numerical 
integrations  only  for  the  symmetry  unique  atoms,  as  we  can  see  from  the  timing  data  the 
symmetry  processing  in  the  current  ACES  II  KS-DFT  gradient  implementation  needs  to  be 
improved  to  take  advantage  of  the  full  savings  offered  by  symmetry.  This  improvement 
requires  computing  the  gradients  in  symmetry  adapted  nuclear  center  basis.  The  efforts 
are  currently  underway  to  improve  upon  the  way  symmetry  being  handled  in  ACES  II  KS- 
DFT  implementation.  However,  the  determination  of  rate  constants  for  larger  molecules 
seldom  can  use  symmetry,  so  for  most  cases  our  experience  with  the  unsymmetrical  N10 
would  be  indicative. 
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Table  8:  A  comparison  of  timing  data  (in  seconds)  of  ACES  II  DFT  gradients  with 
NWChem  (timing  data  are  for  the  SunOs  operating  system) 


C6H6  (cc-pVDZ,  114  AOs) 
B2H6  (cc-pVTZ,  144  AOs) 
N10  (cc-pVDZ,  140  AOs) 


ACES  II 

NWChem 

213 

45 

213 

96 

862 

1011 

It  is  well  known  that  the  accuracy  of  the  subsequent  rate  constants  are  for  the  most  part 
dictated  by  the  quality  of  the  electronic  structure  data.  While  it  is  possible  to  use  the  state 
of  the  art  methods  such  as  CCSD(T)  with  large  basis  sets  to  obtain  highly  accurate  data, 
the  cost  involved  in  such  a  calculation  is  not  practical  in  some  cases.  One  such  case 
would  be  when  the  reactions  of  interest  involve  large  molecules  an  example  of  which 
would  be  reactions  that  are  of  interest  to  biochemists.  Another  such  case  is  when  there 
are  many  reactions  that  the  rate  constants  need  to  be  computed  for  as  in  the  case  of  an 
atmospheric  decomposition  of  a  pollutant  via  a  multitude  of  paths.  ACES  III,  which  is  the 
newly  written,  massively  parallel  version  that  succeeds  ACES  II,  ameliorates  this  partly  by 
permitting  CCSD(T)  and  ACCSD(T)  calculations  to  be  done  on  thousands  of  processors. 


Other  Improvements  to  the  ACES  II  Program  System: 

It  was  observed  that  as  the  system  size  gets  bigger  (especially  for  cluster  and  rings),  the 
geometry  optimization  with  user  defined  internal  or  redundant  internals  becomes 
impractical:  generating  user  defined  internal  input  is  cumbersome  when  large  number  of 
degrees  of  freedoms  are  involved  and  the  automatic  generation  of  internals  done  in  the 
RIC  scheme  does  not  have  adequate  controls  to  limit  the  number  of  internals  that  it 
generates.  In  order  to  remedy  these  problems  ACES  ll/lll  geometry  optimization 
algorithms  were  enhanced  to  work  with  pure  Cartesian  coordinates  (not  recommended  for 
small  symmetric  molecules).  Also,  an  option  to  specify  the  connectivities  and  bypass  the 
automatic  generation  of  internal  coordinates  to  limit  the  number  of  internals  to  the  user’s 
requirements  is  added  to  the  RIC  generation  algorithm.  Another  possible  solution  to  the 
RIC  over  assignment  problem  is  to  prescreen  the  internal  that  it  automatically  assigns  and 
then  eliminate  the  ones  that  are  undesired.  Having  such  an  alternative  would  minimize  the 
user  involvement  and  is  being  currently  implemented. 
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Development  of  ACESRATE  Program  and  Rate  Constant  Calculations 

In  order  to  fully  utilize  the  new  capabilities  being  added  to  ACES  II  and  its  parallel 
version  ACES  III,  we  have  developed  a  program  capable  of  computing  rate  constants 
with  well  established  methods  such  as  transition  state  theory  (TST)  and  variational 
transition  state  theory  (VTST)  in  house. 

CH+N2  Reaction: 


The  reaction  of  CH  with  N2  is  of  considerable  importance  in  high  temperature 
combustion  chemistry27.  We  have  applied  ACES  II  for  the  electronic  state  calculations 
and  ACESRATE  for  the  IRC  and  rate  constant  calculations  of  the  reaction.  We 
proposed  the  new  reaction  mechanism  involving  the  unreported  intermediate,  HNNC, 
and  investigated  the  validity  of  the  mechanism  by  calculating  the  rate  constants  and  by 
comparing  them  to  experiment28.  The  reaction  consists  of  following  nine  elementally 
reactions’ 

(1)  CH  +  N2  — ►  INT 1  (2)  INT 1  — *•  INT2  (3)  INT2  -»  INT3 

(4)  INTI  ->  INT4  (5)  INT4  INT3  (6)  INT2  -►  HNNC 

(7)  HNNC  -» INT4  (8)  INT3  —  NCN  +  H  (9)  HNNC  CNN  +  H 

The  reaction  mechanism  is  shown  schematically  in  figure  5.  The  reactions  (1)  to  (7) 
are  with  corresponding  transition  states  and  the  TST  method  are  used  for  the  rate 
constant  calculation,  while  for  the  reactions  (8)  and  (9)  which  are  dissociations  without 
TS’s,  the  VTST  rate  constants  are  calculated  by  searching  for  the  maximum  of  the  free 
energy  (AG)  potential.  The  geometric  and  electronic  structures  are  calculated  by 
means  of  the  CCSD(T)  level  of  theory  with  the  large  basis  set,  aug-cc-pVTZ. 
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Figure  5.  Schematic  representation  of  the  reaction  mechanism  calculated  at  the 
CCSD(T)/aug-cc-pVTZ  level. 


Table  9  shows  the  calculated  rate  constants  of  the  respective  elemental  processes. 
Within  this  reaction  mechanism,  the  reaction  (8)  can  be  considered  to  be  the  rate 
determining  step.  Although  the  agreement  with  the  experimental  value  is  reasonable 
and  better  at  higher  temperatures,  the  calculated  rate  constant  is  still  smaller  by  about 
an  order  of  magnitude  at  the  highest  temperature.  However,  it  is  known  that  more 
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accurate  electron  correlation  treatments  and  larger  basis  sets  reduce  the  barrier  height 
of  the  reaction,  to  give  a  larger  rate  constant.  We  will  report  basis  set  limit  CCSD(T) 
thermodynamics  and  kinetics  of  CH2+N2  reaction  in  a  review  article  dedicated  to  this 
reaction"9. 


Table  9.  Rate  constants  for  the  respective  elementally  reaction  calculated  by 
CCSD(T)/aug-cc-pVTZ.  The  experimental  values  are  also  shown. 


reaction 

k  [cm3  mol 1  s  ’] 

2662K 

3065K 

3352K 

3550K 

3820K 

(1) 

5.313E+12 

8.455E+12 

1.121E+13 

1.338E+13 

1.668E+13 

(2) 

1.013E+13 

1.142E+13 

1.223E+13 

1.273E+13 

1.337E+13 

(3) 

1.525E+11 

3.881  E+11 

6.590E+1 1 

9.036E+1 1 

1.319E+12 

(4) 

2.723E+08 

1.617E+09 

4.433E+09 

8.089E+09 

1.662E+10 

(5) 

1.455E+13 

1.725E+13 

1.900E+13 

2.012E+13 

2.156E+13 

(6) 

1.239E+10 

3.974E+10 

7.689E+10 

1.140E+11 

1.826E+11 

(7) 

1.570E+10 

4.397E+10 

7.876E+10 

1.115E+11 

1.689E+11 

(8) 

8.491  E+07 

6.674E+08 

2.131  E+09 

4.245E+09 

9.653E+09 

(9) 

4.134E+08 

2.580E+09 

7.252E+09 

1.340E+10 

2.791  E+10 

exptl. 

5.63E+10 

1.30E+11 

1.40E+11 

2.00E+1 1 

2.55E+1 1 

Thermal  Rate  Constants  of  a  Selection  of  Gas  Phase  Reactions  from  the 
Calibration  set  Established  by  Truhlar  and  Coworkers. 

The  electronic  structure  calculations  for  all  40  reactions  in  the  set  are  carried  out  with 
CCSD(T)/MG3,  ACCSD(T)/MG3  and  HF-DFT/MG3  level  to  obtain  reactants  and 
transition  states  structures,  energies  and  the  vibrational  frequencies  in  order  to  obtain 
the  corresponding  TST  and  VTST  rate  constants  (the  frequencies  oat  extra  points  near 
the  transition  state  are  also  obtained  at  the  corresponding  level).  A  subset  from  the 
entire  collection  of  results  is  presented  here  for  the  discussion  purposes  and  the  entire 
collection  will  be  published  elsewhere31,32. 

CCSD(T)  and  ACCSD(T) 

In  Table  10,  we  show  the  TST  rate  constants  obtained  at  the  CCSD(T)  level.  We  also 
calculate  the  full  Arrhenius  plot  and  compare  it  with  the  literature  values.  Figure  6  shows 
the  Arrhenius  plots  for  three  of  the  hydrogen  exchange  reactions  at  the  CCSD(T)  level. 
We  observe  a  generally  good  agreement  in  the  shape  of  the  curve  and  the  values 
between  calculated  rate  constants  and  experimental  rate  constants  for  the  temperature 
range  >  800  K.  Also  note  that  when  more  than  one  experimental  plot  is  available,  the 
theoretical  curves  may  be  used  to  arbitrate  discrepancies  between  sets  of  experimental 
data.  Another  that  at  different  temperatures,  different  transition  states  may  be  accessed 
for  a  particular  reaction  which  is  usually  indicated  by  good  agreement  with  the 
calculated  rate  constant  in  one  temperature  range  and  poor  agreement  at  a  different 
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temperature  range.  Therefore  imperfect  agreement  between  experiment  and  calculation 
may  not  always  be  due  to  the  poor  quality  of  the  electronic  structure  calculation 

Table  10:  Rate  constants  computed  at  700  K.  This  temperature  was  chosen  as 
most  experimental  data  in  the  literature  has  data  points  in  this  temperature 
range 

reactions  E  xp.  Lit.  Value  1  E  xp.  Lit.  Value  2  C  C  S  D(T)/MG3  Abs olute  E  rror 


Cl  +  H2  ->  HCI+  H 

1.38E-012 

1.79E  -012 

1.61E-012 

2.60E  -14 

H  +  HC 1  — >  Cl  +  H2 

2.02E  -012 

1.40E-012 

1  24E-012 

4  70E  13 

OH  +  H2  ->  H  +  H20 

6.28E-013 

5.66E  -013 

1.02E-014 

5  87E-13 

H  +  H20  ->  OH  +  H2 

2.49E-017 

4.71E  *017 

2.96E-017 

6  35E  18 

F  +  H2  ->  H  +  HF 

7.52E  -Oil 

3.71E-010 

2.96E-10 

H  +  HF  *>  F  +  H2 

9.34E-021 

2.96E  -020 

2.03E  20 

OH  +  H  ->  H2  +  0 

4.62E-014 

1.47E  013 

1.08E  014 

8.58E-14 

0  +  H2  ->  0  H  +  H 

3.68E-014 

3  42E-014 

3.81E-015 

3.17E  -14 

Figure  6  (A,  B,  C):  Temperature  dependence  of  rate  constants  for  three 
hydrogen  exchange  reactions  selected  from  Table  3  at  the  CCSD(T)/MG3  level 
of  optimization  and  frequency  calculations. 
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B.  H  +  HF  — >  F  +  H2 


Temperature  Dependence 

H2  +  CI  HCI  +  H 
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C.  OH  +  H  —  H2  +  0 


Temperature  Dependence 


OH  +  H  H2  +  O 


HF-DFT 

A  selection  of  HF-DFT  TST  rate  constants  computed  at  600K  are  shown  in  Table  1 1 . 
For  those  reactions  that  are  shown  in  Table  1 1 .  the  HF-DFT  TST  rate  constants  show  a 
remarkable  agreement  with  experiment.  For  reactions,  F  +  H2  HF,  0  +  H2  O  +  H2, 
OH  +  H  H2  +  O  we  note  that  the  HF-DFT  rate  constant  are  in  agreement  with  the 
CCSD(T)  values  shown  Table  10.  Such  comparison  for  all  the  reactions  in  the  list  will 
be  published31 

Table  11:  HF-DFT  rate  constants. 


Reactions 

Experimental  forward  rate 

HF-DFT  forward  rate 

constant  (cm3  molecule^s  1 

constant  (cm1  molecule  1s'1) 

at  600K 

at  600K 

ch3 

+  h2 

H  + 

1.7*10  lb 

1 .35*1  O’ y 

H 

+  ch4 

CH3+  H \^~ 

2.57x10T3 

6.26*1 0715- 

H 

+  ph3 

h2  +ph2 

6.2*10 

1 .15*10  11 

F 

+  h2 

HF  +  H34- 

7.39x10  11 

2.52*10  11 

0 

+  h2 

OH  +H54- 

1 .14x10  44 

4.39*10^ 

OH 

+  H 

0  +  H^4 

1.89x10  14 

1 .37*1 0"^- 

0 

+  HCI 

OH  +  Cl 34 

2.62x10  T4 

2.68*10  15 

l0H 

+  Cl 

0  +  hcF4- 

8.34x1 044- 

3.00*1 0747"- 

Summary. 

New  methods  including  HF-DFT  and  ACCSD(T)  have  been  developed  in  this  STTR, 
including  their  analytical  gradients.  We  have  also  developed  a  much  improved  method 
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for  generating  accurate  IRC’s.  In  addition  we  have  interfaced  ACES  II  and  III  to  our  own 
ACESRATE,  also  developed  for  this  project,  and  other  packages  like  POLYRATE  to  do 
such  rate  constant  calculations  while  using  the  electronic-structure  information  from 
ACES.  The  latter’s  performance  was  enhanced  while  adding  analytical  gradients  for 
DFT  and  the  previously  unexplored  HF-DFT,  which  is  the  simplest  possible  correlated 
method  for  determining  the  electronic  structure  of  molecules.  This  is  complimented  by 
ACCSD(T)  which  is  the  highest  level  of  electronic  structure  theory  that  can  be  expected 
to  be  routinely  applicable  to  a  large  number  of  systems.  Extensive  numerical  results 
show  the  performance  of  these  tools  for  rate  constants  and  related  properties.  As  ACES 
II,  ACES  III,  and  ACESRATE  are  freely  distributed,  these  methods  can  be  extensively 
applied  by  non-experts  to  their  own  problems. 
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